# 计算 E(|X-a|) 其中X服从正态分布

.exp = function(mu, sigma, a)
    return(exp(-(a - mu) ^ 2 / (2 * sigma ^ 2)) *
               sqrt(2 / pi) * sigma + (a - mu) *
               VGAM::erf((a - mu) / (sigma * sqrt(2))))

test = function (mu, sigma, a, n = 1e6) {
    dat = rnorm(n, mu, sigma)
    cat(paste0(
        "simulation: E(|X",
        sprintf("%+.6f", a),
        "|)= ",
        sprintf("%.6f", mean(abs(dat - a))),
        "\n"
    ))
    cat(paste0(
        "     exact: E(|X",
        sprintf("%+.6f", a),
        "|)= ",
        sprintf("%.6f", .exp(mu, sigma, a)),
        "\n"
    ))
}

for (i in 1:10) {
    mu = runif(1, -10, 10)
    sigma = runif(1, 0, 10)
    a = runif(1, -10, 10)
    n = 4e7
    cat(sprintf("X ~ N(%.6f, %.6f)\n", mu, sigma))
    test(mu, sigma, a, n)
    cat("\n")
}
